The phase diagram of numerous materials of technological importance features high-symmetry high-temperature phases that exhibit phonon instabilities. Leading examples include shape-memory alloys, as well as ferroelectric, refractory, and structural materials. The thermodynamics of these phases have proven challenging to handle by atomistic computational thermodynamic techniques, due to the occurrence of constant anharmonicity-driven hopping between local low-symmetry distortions, while maintaining a high-symmetry time-averaged structure. To compute the free energy in such phases, we propose to explore the systems potential-energy surface by discrete sampling of local minima by means of a lattice gas Monte Carlo approach and by continuous sampling by means of a lattice dynamics approach in the vicinity of each local minimum. Given the proximity of the local minima, it is necessary to carefully partition phase space by using a Voronoi tessellation to constrain the domain of integration of the partition function, in order to avoid double-counting artifacts and enable an accurate harmonic treatment near each local minima. We consider the bcc phase of titanium as a prototypical examples to illustrate our approach.