Summary A Bayesian statistical framework was previously developed for modal identification of well-separated modes incorporating ambient vibration data, that is, operational modal analysis, from multiple setups. An efficient strategy was developed for evaluating the most probable value of the modal parameters using an iterative procedure. As a sequel to the development, this paper investigates the posterior uncertainty of the modal parameters in terms of their covariance matrix, which is mathematically equal to the inverse of the Hessian of the negative log-likelihood function evaluated at the most probable value. Computational issues arising from the norm constraint of the global mode shape are addressed. Analytical expressions are derived for the Hessian so that it can be evaluated accurately and efficiently without resorting to finite difference. The proposed method is verified using synthetic and laboratory data. It is also applied to field test data, which reveals some challenges in operational modal analysis incorporating multiple setups. Copyright © 2014 John Wiley & Sons, Ltd.