  We develop a continuous-time quantum Monte Carlo method based on a
strong-coupling expansion for Anderson impurity models with phonon-assisted
hybridizations for arbitrary number of phonon modes. As a benchmark, we
investigate the two-channel Anderson model with a single phonon, and
numerically demonstrate that an SO(5) susceptibility composed of
localized-electron charge and phonon-parity operators diverges logarithmically
at the non-Fermi liquid critical point in the model, which verifies the
predictions by the boundary conformal field theory[K. Hattori: Phys. Rev. B
{\bf 85} (2012) 214411].
