Monte Carlo algorithms are barely considered in spin foam quantum gravity. Due to the quantum nature of spin foam amplitudes one cannot readily apply them, and the present sign problem is a threat to convergence and thus efficiency. Yet, ultimately the severity of the sign problem in spin foams is not known. In this article we propose a new probability distribution for coherent (boundary) intertwiners, which we use to define a Markov Chain Monte Carlo algorithm. We apply this algorithm to the SU(2) coherent vertex amplitude for various Regge-type boundary data and find convergent, accurate results at far lower costs than the explicit calculation. The resources are instead used to increase the size of boundary spins, bridging the gap to the asymptotic formulae. While the sign problem is not solved, it is under control in the vast majority of cases. We close by discussing how this algorithm can be extended to larger triangulations with boundary and 4d Lorentzian spin foam models. More speculatively, we discuss how this algorithm can be used to sample (Regge-geometric) bulk intertwiners.

42 pages, 15 figures; v2: matching version accepted for publication in PRD, expanded on explanation of algorithm, added appendix describing convergence properties of algorithm for different boundary data