We introduce variational wave functions to evaluate the ground-state properties of spin-phonon coupled systems described by the Su-Schrieffer-Heeger model. Quantum spins and phonons are treated on equal footing within a Monte Carlo sampling, and different regimes are investigated. We show that the proposed variational Ansatz yields good agreement with previous density-matrix renormalization group results in one dimension and is able to accurately describe the spin-Peierls transition. This variational approach is constrained neither by the magnetoelastic-coupling strength nor by the dimensionality of the systems considered, thus allowing future investigations in more general cases, which are relevant to spin-liquid and topological phases in two spatial dimensions.
