We revisit the problem of inferring mantle viscosity from postglacial relative sea level (RSL) data across the Hudson Bay. We invert a recently revised data set using the Metropolis algorithm together with an annealing schedule: this method, which is well established in geophysics, is applied here for the first time to the glacial isostatic adjustment problem. The Metropolis algorithm performs a search, which is not limited to downhill moves in the model space and thus is less influenced by local minima of the misfit than traditional inverse approaches. Furthermore, its CPU requirements are far superior to Monte Carlo methods. The major drawbacks include slow convergence and the need for careful tuning of crucial variables such as the temperature schedule and the increment in the model space. When all the Hudson Bay RSL data are considered, and the viscosity of the upper mantle above the 670 km discontinuity is inverted, the best fitting solution is characterized by a viscosity close to 2 ¿ 1020 Pa s. However, when the shallow upper mantle and transition zone viscosity are separately inverted, other less traditional solutions with a more complex viscosity structure are found to be equally possible. A stable feature is the lower mantle viscosity, which is generally found to be close to the value of 1021 Pa s in all of the stochastic inversions we have performed. The solutions agree with previous findings concerning both postglacial rebound observables and global geodynamics signatures.