  By combining classical Monte Carlo and Bethe ansatz techniques we devise a
numerical method to construct the Truncated Generalized Gibbs Ensemble (TGGE)
for the spin-1/2 isotropic Heisenberg ($XXX$) chain. The key idea is to sample
the Hilbert space of the model with the appropriate GGE probability measure.
The method can be extended to other integrable systems, such as the
Lieb-Liniger model. We benchmark the approach focusing on GGE expectation
values of several local observables. As finite-size effects decay exponentially
with system size, moderately large chains are sufficient to extract
thermodynamic quantities. The Monte Carlo results are in agreement with both
the Thermodynamic Bethe Ansatz (TBA) and the Quantum Transfer Matrix approach
(QTM). Remarkably, it is possible to extract in a simple way the steady-state
Bethe-Gaudin-Takahashi (BGT) roots distributions, which encode complete
information about the GGE expectation values in the thermodynamic limit.
Finally, it is straightforward to simulate extensions of the GGE, in which,
besides the local integral of motion (local charges), one includes arbitrary
functions of the BGT roots. As an example, we include in the GGE the first
non-trivial quasi-local integral of motion.
