  We derive a semi-analytical form for the Wigner transform for the canonical
density operator of a discrete system coupled to a harmonic bath based on the
path integral expansion of the Boltzmann factor. The introduction of this
simple and controllable approach allows for the exact rendering of the
canonical distribution and permits systematic convergence of static properties
with respect to the number of path integral steps. In additions, the
expressions derived here provide an exact and facile interface with quasi- and
semi-classical dynamical methods, which enables the direct calculation of
equilibrium time correlation functions within a wide array of approaches. We
demonstrate that the present method represents a practical path for the
calculation of thermodynamic data for the spin-boson and related systems. We
illustrate the power of the present approach by detailing the improvement of
the quality of Ehrenfest theory for the correlation function
$\mathcal{C}_{zz}(t) = \mathrm{Re}\langle \sigma_z(0)\sigma_z(t)\rangle$ for
the spin-boson model with systematic convergence to the exact sampling
function. Importantly, the numerically exact nature of the scheme presented
here and its compatibility with semiclassical methods allows for the systematic
testing of commonly used approximations for the Wigner-transformed canonical
density.
